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Abstract. One can solve the Jeans equation analytically for equilibrated dark matter 
structures, once given two pieces of input from numerical simulations. These inputs are 
1) a connection between phase-space density and radius, and 2) a connection between 
velocity anisotropy and density slope, the a — (3 relation. The first (phase-space 
density v.s. radius) has already been analysed through several different simulations, 
however the second {a — /3 relation) has not been quantified yet. We perform a large 
set of numerical experiments in order to quantify the slope and zero-point of the 
a — f3 relation. We find strong indication that the relation is indeed an attractor. 
When combined with the assumption of phase-space being a power-law in radius, this 
allows us to conclude that equilibrated dark matter structures indeed have zero central 
velocity anisotropy (3q « 0, central density slope of ~ —0.8, and outer anisotropy 
of /3oo « 0.5. 



The a — (3 relation 



2 



1. Introduction 

We have seen remarkable progress in the understanding of pure dark matter structures 
over the last few years. This was triggered by numerical simulations which have 
observed general trends in the behaviour of the radial density profile of equilibrated 
dark matter structures from cosmological simulations, which roughly follow an NFW 
profile P la El il 13 El, see also 13 GDI HB 112 ESI for references. General trends 
in the radial dependence of the velocity anisotropy have also been suggested [HI ITKj . 
Recently, more complex relations have been identified, holding even for systems that do 
not follow the simplest radial power-law behaviour in density. These relations are first 
that the phase-space density, p/cr^, is a power-law in radius jTH], and second that there 
is a linear relationship between the density slope and the anisotropy jTTj. A connection 
between the shape of the velocity distribution function and the density slope has also 
been suggested [13 HH] ■ 

Using the Jeans equation together with the fact that phase-space density is a power- 
law in radius allows one to find the density slope in the central region numerically |16j . 
and even analytically for power-law densities jSO]. These results were generalized in |21j . 
Recently a more refined study [22] using both the phase-space density being a power- law 
in radius and also the linear relationship between density slope and anisotropy, showed 
that one can solve the Jeans equations analytically and extract the radial dependence 
of density, anisotropy, mass etc. 

It therefore appears that we only have to understand the two numerical relationships 
described above in order to fully quantify pure dark matter halo structures. However, 
no theoretical explanation for the origin of these relations is known to us. 

The relationship between phase-space density and radius has been considered 
several times [13 123 1211 1221 123 1211 12S1 and seems to be well established, however, the 
other crucial ingredient in the analysis, namely the linear relationship between density 
slope and anisotropy has only been investigated qualitatively [T71. In this paper we 
perform a large set of simulations in order to quantify this relationship. We show that 
with present day simulations there does indeed appear to be a linear relation between 
density slope and anisotropy. When combined with the assumption of phase-space being 
a power-law in radius this implies zero anisotropy near the center with density slope of 
approximately —0.8, and that the outer anisotropy is radial and close to +0.5. 

2. Head-on collisions 

The first controlled simulation is the head-on collision between two initially isotropic 
NFW structures. For the construction of the initial structures, we use the Eddington 
inversion method as described in |26j. We create an initial NFW structure with zero 
anisotropy containing 1 million particles. The parameters are chosen to correspond to 
a total mass of 10^^ solar masses (concentration of 10), with half of the particles within 
160 kpc. The structure is in equilibrium in the sense that when evolving such a structure 
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in isolation its global properties remain unchanged. 

We now place two such structures very far apart with 2000 kpc between the centres, 
which is well beyond the virial radius. Using a softening of 0.2 kpc, we now let these 
two structures collide head-on with an initial relative velocity of 100 km/sec towards 
each. After several crossings the resulting blob relaxes into a prolate structure. We 
run all simulations until there is no further time variation in the radial dependence of 
the anisotropy and density. We check that there is no (local) rotation. We run this 
simulation for 150 Gyrs (corresponding rougly to 10 Hubble times), which means that 
a very large part of the resulting structure is fully equilibrated. 

All simulations were carried out using PKDGRAV, a multi-stepping, parallel 
code |2Z1, which uses spline kernel softening and multi-stepping based on the local 
acceleration of particles. Force accuracy is set by an opening angle of B = 0.7 which 
combined with the use of a 4-th order multipole expansion results in typical RMS force 
errors of better than 0.1%. The simulations were performed on the zBox and zBox2 at 
the University of Zurich %. 

We now extract all the relevant parameters in radial bins, logarithmically 
distributed from the softening length to beyond the region which is fully equilibrated. 
The resulting density profile is very similar to an NFW profile. We calculate the radial 
derivative of the density (the density slope) 

and the velocity anisotropy 

(2) 

r 

in each radial bin, where the a^j, are the one dimensional velocity dispersions in the 
tangential and radial directions respectively. This is shown in figure ^ as green filled 
pentagons. We also present the initial conditions as a red (straight) line. It is clear from 
this figure that some connection between a and (3 exists, in the sense that a density 
slope of about -1 has small anisotropy, and density slope of about -3 has large radial 
anisotropy. In these figures we are including points which are inside the resolved region 
and also points which are outside the fully equilibrated region. We will discuss this issue 
fully in section 01 

The first important question to address is if the a — (3 relation found after this 
collision represents an optimal configuration of the dark matter system, or if subsequent 
collisions might lead to a different configuration. We therefore take the resulting 
structure (now containing 2 million particles, in a prolate shape) and collide 2 copies 
together, again separated by 2000 kpc, and with initial relative velocity of 100 km/sec. 
We now use softening of 1 kpc to speed up the simulation. The resulting structure is 
even more prolate when observed in density contours, and we evolve this for another 
150 Gyr. It should be noted that already after a few crossing times the structures are 

I |http : //www-theorie -physik.unizh. ch/~stadel/zBo3{| 
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Figure 1. Head-on collision between two initially isotropic NFW structures. The red 
(straight) line shows the initial condition for each of the colliding structures. The green 
filled pentagons show the result after the first collision, where the structure contains 
2 million particles. The blue open triangles show the result after a second collision (4 
million particles) . The green and blue symbols land very near each other in the region 
inside of a density slope of about —2.2 



sitting near the a — (3 relation, and the fact that we run the simulations for much longer 
times has only a small effect in the outer region, where the particle density is much 
smaller. 

We again extract in radial bins, and the a — (3 is shown as blue open triangles in 
figure^ We see very clearly, that for density slopes more shallow than about —2.2, there 
is virtually no change in the a — j3 relation. We conclude that in the central region, 
where a significant perturbation occurred during the collision, the a — (3 relation is 
unchanged, and hence the structure has indeed reached an optimal state. In the outer 
region there has been a small change in the connection between a and /?, in the sense 
that the first collision brought the outer region away from the initial conditions and the 
second collision brought the system even further away. 
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2.1. Dependence on shape? 

The resulting structures described above are prolate when observed in density contours. 
One may fear that the resulting a — j3 relation will depend strongly on the axis-ratios 
of these structures. Naturally, if the relevant axis-ratios are instead the ones extracted 
when considering contours in potential energy, then we should expect only an effect in 
the very central region, since the potential energy contours are close to spherical in the 
outer region. 
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Figure 2. Repeated head-on collision between two initially isotropic NFW structures, 
to test the effect of shape. The different colours refer to the axis along which the 
second collision was made. In the inner region there is virtually no difference in the 
resulting a — (3 relations. 



In order to test this question, we extracted the resulting prolate structure after the 
first collision described above. We can now decide to collide this structure along different 
axes, either along the long, intermediate or short axes. We perform 2 test collisions, 
one is a second collision along the intermediate axis, and the other is a second collision 
along the short axis (there is very little difference between the short and intermediate 
axis after the first collision). These structures end up being triaxial, almost oblate. 
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The softening of these last coUisions was chosen to be 2 kpc in order to speed up the 
simulation. The collision along the long axis was described above, resulting in a very 
prolate structure. The results are shown in figure El 

It is clear from figure El that in the inner region (inside density slope of about 
—2.2) there is very little difference between the 3 different structures. These 3 resulting 
structures have rather different shape when seen in density contours, but are very 
similar when seen in potential energy contours. We conclude from this that whereas the 
definition of density slope does depend slightly on the shape of the density contours, 
then the a — (3 relation is almost independent of the shape j2H] • 

We note that the structures should be more perturbed when collided along the long 
axis, and we do indeed observe that the changes in the a — (3 relation in the outer region 
(beyond density slope of —2.2) are larger when collided along the long compared 
to collisions along the short axis. 

2.2. Dependence on initial conditions? 

To address the question of how strongly the a — (3 relation from the head-on collisions 
described above depend on the initial conditions, we now consider 2 simulations each 
with two steps, la) First we create a spherical isotropic NFW structure with 1 million 
particles as described above. We take each individual particle in this structure and 
put its total velocity along the radial direction, maintaining the sign with respect to 
inwards or outwards from the halo centre. We keep the energy of each particle fixed. 
This structure is thus strongly radially anisotropic, lb) Now we take two such radially 
anisotropic structures and place their centres 2000 kpc apart (well outside the virial 
radius) with relative velocity of 100 km/sec towards each other. We use softening of 1 
kpc for these tests. Before the centres of the two structures collide the radial motion 
of the particles behaves almost like a radial infall simulation, and strong scattering 
of the particles is observed. The 2 individual structures pick random orientations in 
space as they become triaxial j^; their orientations are completely erased after the 
two structures collide head-on. We let the structure equilibrate. We take two copies 
of the resulting structure and collide these head-on, again with 100 km/sec. The final 
equilibrated region contains approximately 2.5 million particles (of the total of 4 million 
particles). 

We follow this with a simulation identical to the one just described, except 
that we place each individual particle in the initial structures on a tangential orbit, 
without changing the energy of the individual particles. 2a) This initial condition thus 
corresponds to strongly tangential motion. When letting this structure equilibrate in 
isolation it settles down with a tangential anisotropy of /3 ~ —2. 2b) Again, two copies 
of this equilibrated structure are collided head-on as in step lb. The result of these tests 
are shown in figure El 

We see from figure El that the resulting a — 13 relation is virtually independent of 
initial conditions, and we conclude that the collisions were sufficiently violent to erase 
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Figure 3. Repeated head-on collision between two NFW structures, to test the effect 
of different initial conditions. The different colours refer to different initial conditions; 
green (open) circles show the results from the initially radially anisotropic case (1), 
cyan (filled) pentagons show the results of the initially tangentially anisotropic case 
(2), and blue (open) triangles was the isotropic collision described in section |21 In the 
region inside of a density slope of about —2.4, there is virtually no different between 
the resulting a — (3 relations. 



the initial conditions sufficiently to conform with the a — (3 relation. We emphasize, 
that we are not stating that initial conditions are erased completely, however, that they 
are erased sufficiently to allow the resulting structure to obey the a — (3 relation. 

2.3. Symmetric simulations 

The simulations described above were all performed through very non-symmetric 
collisions. We therefore performed two experiments with symmetric collisions. First 
we take 6 copies of the isotropic spherical NFW structure, each with 1 million particles, 
and place them symmetrically along the x, y and z axes. We collide these with 
initial separation and relative velocities similar to what was described above. For this 
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simulation we use softening of 2 kpc. The resulting structure has 6 million particles. 
Second, we created an NFW structure with the same parameters as described above, 
but only containing 10^ particles. We now place 15 of these in a cubic symmetry (1 
cell-centered, 6 face-centered and 8 at corners) with initial separation about 600 kpc 
and initial relative velocities of 60 km/sec, and collide these together using softening of 
1 kpc. After letting the resulting structure equilibrate we take this resulting structure, 
place 15 copies in a cubic symmetry, and collide these again. The resulting structure 
thus contains 2.25 million particles. The results of these two simulations are shown in 
figure m There is good agreement between the resulting a — (3 relation for these two 
structures in the central region. 
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Figure 4. Symmetric collision between 6 large isotropic NFW structures (green open 
triangles); repeated collisions of numerous small initially isotropic structures (red filled 
pentagons); tangential instability (black open circles). These significantly different 
simulations exhibit remarkably similar a — p relations. 
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2.4- Tangential instability 

One may fear that the resuhing structures from the various simulations shown thus far 
end up near the same hne in a — f3 space simply because they sample only relatively 
similar kinds of merging processes. The concern is that all the global properties of the 
final structures are the result of the fact that we slam structures into one another. Also 
the radial infall simulations included in ref. ^Tj can be argued to be nothing but small 
structures falling into some larger structure. 

In order to probe this issue further, we will simulate the evolution of a rather 
artificial structure. We populate a small number of particles, 4 ■ 10^, uniformly (but 
with randomly chosen positions) on an infinitely thin spherical shell of radius 0.01 in 
simulation units. These particles all have zero initial velocity. Then we place 1 million 
particles on another infinitely thin spherical shell at the radius 20 in simulation units, 
again uniformly distributed, but with randomly chosen positions. These particles are 
all placed on exactly circular orbits, but with random directions (tangentially to this 
sphere). All particles have exactly the same mass. Now, if this system had infinitely 
many particles uniformly distributed then the system would be in equilibrium, however, 
the poissonian noise of the particles in the initial condition, and the numerical noise 
from the integration of the equation of motion, imply that the particles will slowly 
start clumping together at random points on the sphere. The central collection of 
particles (4 ■ 10^), which effectively act as a central massive mass point, quickly settles 
into a small equilibrated structure, leaving the gravitational influence on the outer 
particles virtually unaffected. The particles on the outer shell, however, slowly break 
the initial spherical symmetry. Particles passing nearby each other start scattering off 
each other, kicking some particles outwards and some particles inwards. Eventually the 
entire system breaks up, and the small initial central collection of particles are swallowed 
by the large collection of particles from the outer shell. 

We run this simulation with softening of 0.1 in these simulation units, and we 
let the system evolve towards a new equilibrium state. This state does, surprisingly 
enough, resemble an NFW structure in density, and we extract the a — (3 relation for 
this resulting structure. The results are plotted as black (open) circles on figure EJ and 
the agreement with the other simulations is striking. 

We conclude that the appearance of the a — P relation must have a deeper 
foundation than merely being the result of considering systems which were constructed 
to smash into each other. 

2.5. Comparison with previous works 

In order to compare the results of these various simulations presented above, we now 
plot the 3 double head-on collisions with different initial conditions, the two symmetric 
collisions (between 6 and 225 initially isotropic NFW structures), and the tangential 
instability in figure El On the same figure we plot the line which was argued to be 
a reasonable fit in ref. J7]. This line was a fit to simulations including mergers of 
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disk galaxies, cosmological structures, and radial infall simulations, all simulations very 
different from the ones considered in this paper. We see indeed that our simulations are 
in fair agreement with the results of ref. jUj. 




density slope 



Figure 5. A collection of the various simulations described in the previous sections. 
The straight (green) line is from ref. |17j , which was shown to provide a reasonable fit 
to a collection of different simulations, including mergers of disk galaxies, cosmological 
structures, and radial infall simulations. 



3. Region of trust 

In all the figures above we have plotted a very extended region in the density slope, 
which corresponds to plotting points very near the central region (small negative a) 
and very distant (large negative a). 

In the central region one can generally only trust a region beyond a few times the 
softening length. In almost all of our simulations this corresponds to a region where 
the density slope is about -1. Further inwards numerical noise leads to density profiles 
more shallow than -1, which cannot be trusted. However, most of our simulations have 
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been run for very long times (much longer than a corresponding typical cosmological 
simulation) and the numerical noise therefore has a slightly larger reach. We are being 
rather conservative, and ignore the regions inside 5 times the softening. Concerning the 
outer region, we have visually inspected all the simulations, and it is clear that regions 
with slopes more shallow than —3 are fully equilibrated. Looking at figure El one also 
sees that the scatter m. a — (3 becomes large beyond a of —3. We therefore trust the 
region with slopes more shallow than a = —3. This region (excluding points inside 5 
times the softening length) is shown for each simulation in figure El together with a line 
from the equation 



which provides a reasonable fit within this region. We see that the scatter in (3 is only 
about 0.05. 

4. Comparison with theory 

In a recent paper, Dehnen & McLaughlin (2005) showed, that under the two assumptions 
that phase-space density is a power-law in radius, and that there is a linear a — jS 
relation, that the central slope of dark matter structures must be ao = —(7 -|- 10/3o)/9, 
where (3o is the central anisotropy. In other words, the central values of the density 
slope and the anisotropy must be somewhere on the thin (blue, solid) line in figure 
however, from theory alone there is no telling where. 

We can now compare this theoretical result with our numerical findings. Our 
numerical results are approximated by the fat (red, dashed) line. We see that the 
two lines cross near (3 = and a = —7/9, showing that the central part of dark matter 
structures is isotropic, /?o ~ 0, and that the central density slope is indeed ao ~ —7/9. 

Theory also shows that the outer density slope is about 7oo = —31/9 ~ —3.44, 
which when compared with our findings result in (3oo ~ 0.53. Thus, from theoretical 
works we know the innermost and outermost density slopes, and combined with the 
results of this paper we now also know both innermost and outermost anisotropy of 
pure dark matter structures. 

Ref. [221 noticed that when assuming that phase-space density is a power-law in 
radius, then the Jeans equations have a particularly simple form if and only if there is 
a linear a — (3 relation as originally suggested in ref. • Our numerical results indeed 
confirm such linearity, and hence support the use of this particularly simple version of 
the Jeans equation for theoretical studies. 

4.1. Comparison with cosmological simulations 

We clarified in section |21 that the a — (3 relation only holds for systems which have been 
perturbed sufficiently and subsequently allowed to relax. One can imagine setting up 
two equilibrium systems each with zero anisotropy, and colliding these with large impact 
parameter, in which case the system is not perturbed sufficiently, and the resulting 
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Figure 6. Trusted region. We liave excluded points inside a central region 
corresponding to 5 times the softening. We do not trust points ftirther out than a ~ —3. 
The red dashed line is a fit to the points in this region, given by /3 « —0.2 (a + 0.8). 
The thin blue (solid) line shows the theoretical results of ref. for the connection 
between the central density profile and the central anisotropy: from theory alone these 
central values can be anywhere on this solid line. The crossing of the two straight lines 
thus shows the actual values for the central density slope and anisotropy. 



system should not be expected to land near the a — (3 relation. As was shown in 
section even with zero impact parameter, only the densest part of the system (inside 
density slope of about —2.2) will reach the a — (3 relation after only one collision. 
One should keep in mind that the problem here is that setting up a system with zero 
anisotropy is very artificial, and not in agreement with structures formed in cosmological 
simulations. 

Similarly, only the central part of structures formed in cosmological simulations 
should be expected to land on the a — [3 relation. This is because most of the infalling 
matter has never been near the center of the structure and hence has not been perturbed: 
most particles outside half of the virial radius have not been through the center and 
have therefore not been mixed sufficiently. This expectation is indeed confirmed by 



The a — (3 relation 



13 



cosmological simulations (private communications with J. Diemand and C. Power, see 
also the large scatter in figure 3 of ref. [22]) which show that the central regions of 
dark matter structure do land on the a — j3 relation, whereas the outermost regions do 
not quite. The a — j3 relation has also been shown to hold for high-cr subset of galaxy 
haloes [HH] . 

5. Attractor solution? 

It appears from the agreement between the different simulations presented above, that 
the a — (3 relation is some kind of attractor. In order to test such a claim one should 
perform minor perturbations of the system, and then observe in which direction the 
system fiows. Antonov's laws of stability tell us that many systems are in equilibrium, 
e.g. an isotropic Hernquist structure will remain isotropic when exposed to minor 
perturbations. Antonov's laws of stability are valid under the assumption that there 
is a zero on the r.h.s. of the Boltzmann equation, and we will therefore expose our 
system to perturbations which act like an instantaneous non-zero term on the r.h.s. of 
the Boltzmann equation. 

We take an isotropic NFW structure which is in total equilibrium. Then we take 
each individual particle and change its velocity with a random number, keeping its 
direction unchanged. This is done in such a way that the new energy is changed by 
maximum 25% from the original energy, and in such a way that there is overall energy 
conservation. This perturbation is thus completely isotropic, and should not by itself 
change the isotropy of the system. Subsequently we let the system relax to a new 
equilibrium state, and we extract density profile and anisotropy of this new system. 

We then take the new equilibrated system and perturb it in a similar way with 
a new set of random numbers. After letting it relax, we can repeat the process. The 
resulting a — /5's are presented in figure [71 where we show the relaxed systems, after 
having perturbed 1, 4, 7 and 10 times. We see indeed, that the system is slowly, but 
surely, moving in the direction of the universal a — (5 relation. We therefore conclude 
that the a — (5 relation really is an attractor. 

We did not perform further perturbations of the system, because for anisotropic 
systems (like after the lO'th perturbation) these isotropic perturbations will on their own 
tend to isotropize the system, and hence the system may potentially appear to not fiow 
towards the attractor. Most probably one must devise more sophisticated perturbations 
in order to see the fiow all the way toward the universal a — (3 line. 

6. Conclusions 

We have quantified the relationship between the density slope and the velocity 
anisotropy, the a — (3 relation. We have performed a large set of simulations to 
investigate systematic effects related to shape and initial conditions, and we find that 
the a — (3 relation is almost blind to both shape and initial conditions, as long as the 
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Figure 7. We perturb the energies sliglitly of each individual particle in an equilibrum 
system. After letting the system relax do we perturb again. We show the resulting 
a — /3's after 1, 4, 7 and 10 perturbations. These perturbations are isotropic, and do 
therefore by themselves not induce any anisotropy. The system does indeed appear to 
move towards the universal a — p relation, indicating that it is indeed an attractor. 



system has been perturbed sufficiently and subsequently allowed to relax. We find strong 
indications that the relation is indeed an attractor. 

We have performed symmetric and highly non-symmetric simulations, along with 
a simulation of the tangential orbit instability to extract the zero point and slope of the 
a — f] relation in the regions which are numerically trustworthy. These simulations 
complement the simulations performed in ref. [12], which included a cosmological 
simulation and radial infall collapse, and yet exhibit a striking level of agreement with 
this work. When compared with analytical results we find that the central region is 
indeed isotropic, and that the outer asymptotic anisotropy is radial, with a magnitude 
of /3Ri0.5. 
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